Atomistic theory for the damping of vibrational modes in mono-atomic gold chains 
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We develop a computational method for evaluating the damping of vibrational modes in mono- 
atomic metallic chains suspended between bulk crystals under external strain. The damping is due 
to the coupling between the chain and contact modes and the phonons in the bulk substrates. The 
geometry of the atoms forming the contact is taken into account. The dynamical matrix is computed 
with density functional theory in the atomic chain and the contacts using finite atomic displacements, 
while an empirical method is employed for the bulk substrate. As a specific example, we present 
results for the experimentally realized case of gold chains in two different crystallographic directions. 
The range of the computed damping rates confirm the estimates obtained by fits to experimental 
data [Prederiksen et at, Phys. Rev. B 75, 205413(R)(2007)]. Our method indicates that an order- 
of-magnitude variation in the harmonic damping is possible even for relatively small changes in the 
strain. Such detailed insight is necessary for a quantitative analysis of damping in metallic atomic 
chains, and in explaining the rich phenomenology seen in the experiments. 

PACS numbers: 63.22. Gh,68.65.-k,73.40. Jn 



I. INTRODUCTION 

The continuing shrinking of electronic devices and the 
concomitant great interest in molecular electronics^ have 
underlined the urgency of a detailed understanding of 
transport of electrons through molecular-scale contacts. 
A particularly important issue concerns the energy ex- 
change between the charge carriers and the molecular 
contact. Thus, the local Joule heating resulting from the 
current passing through the contact, and its implications 
to the structural stability of such contacts are presently 
under intense investigation^^^. Experimentally, local 
heating in molecular conductors in the presence of the 
current has been inferred using two-level fluctuations^ 
and Raman spectroscopy 8 -. 

Mono-atomic chains of metal atoms^ are 
among the simplest possible atomic-scale con- 
ductors. The atomic gold chain is probably the 
best studied atomic-sized conductor, and a great 
deal of detailed information is available from 
experiments^^^^^^^^' 2 ^, and related 
theoretical studiesi^^^^^^ 2 ^^^' 3 ^ 3 .. The 
current induced vibrational excitation and the stability 
of atomic metallic chains have been addressed in a few 
experiments 3 ^^^. 

In the case of a gold chain Agrait et at— reported well- 
defined inelastic signals in the current- voltage character- 
istics. These signals were seen as a sharp 1% drop of the 
conductance at the on-set of back-scattering due to vibra- 
tional excitation when the voltage equals the vibrational 
energy. Especially for the longer chains (6-7 atoms), the 
vibrational signal due to the Alternating Bond-Length 
(ABL) mod o 28 ' 31 , dominates. This resembles the situa- 
tion of an infinite chain with a half-filled electronic band 
where only the zone-boundary phonon can back-scatter 



electrons^ due to momentum conservation. 

The inelastic signal gives a direct insight into how the 
frequency of the ABL-mode depends on the strain of the 
atomic chain. This frequency can also be used to infer 
the bond strength. The signature of heating of the vi- 
brational mode is the non-zero slope of the conductance 
versus voltage beyond the on-set of excitation: with no 
heating the curve would be flat. Fits to the experiment 
on gold chains using a simple model 3 -^ suggest that the 
damping of the excitation, as expected, can be signifi- 
cant. However, the experiments in general show a vari- 
ety of behaviors and it is not easy to infer the extent of 
localization of the ABL vibration or its damping in these 
systems^ 8 -. 

In order to address the steady-state effective temper- 
ature of the biased atomic gold chain theoretically, it is 
necessary to consider the various damping mechanisms 
affecting the localized vibrations, such as their coupling 
to the vibrations in the contact, or to the phonons in 
the surrounding bulk reservoirs. This is the purpose of 
the present paper: we calculate the vibrational modes 
in atomic gold chains and their coupling and the re- 
sulting damping due to the phonon system in the leads. 
We work within the harmonic approximation and employ 
first principles density functional theory (DFT) for the 
atomic chain and the contacts^ while a potential model 
is used for the force constants of the leads 3 -^. 

Experimental TEM studie a 20 i 4Q have shown that 
atomic chains form in the (100) and (111) directions while 
the (110) direction gives rise to thicker rods^. Therefore 
we focus on chains between two (lOO)-surfaces or (111)- 
surfaces. We consider chain-lengths of 3-7 atoms and 
study the behavior of their vibrations and damping when 
the chains are stretched. The TEM micrographs also in- 
dicate that the chains are suspended between pyramids, 
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FIG. 1: Two ways of partitioning the central part of the chain- 
substrate system. [Top] The Chain is the part of the system 
that only contains one atom in a plane parallel to the surface, 
and the Base is what connects the two-dimensional surface to 
the Chain. [Bottom] The Pyramid is the Base plus the Chain 
atom closest to the Base. The Central Chain is the remaining 
part of the Chain after removing one atom at each end. 



so in our calculations we add the smallest possible fcc- 
stacked pyramid to link the chain to the given surfaces. 

As we shall show below, at low strain the gold chains 
have harmonically undamped ABL-modes with frequen- 
cies outside the bulk band. The long chains of 6-7 atoms 
also have ABL-modes with very low damping at high 
strain . Our results indicate that chains between (111)- 
surfaces will have a lower damping than chains between 
(lOO)-surfaces. Importantly, we find that the damping 
is an extremely sensitive function of the external strain: 
an order of magnitude change may result from minute 
changes in the strain. This may provide a key for under- 
standing the rich behavior found in experiments. 

The paper is organized as follows. In Sect. II we 
describe how the central quantities, i.e., the dynamical 
matrix, the projected density-of-states, and the damping 
rates are calculated. Section III is devoted to the analy- 
sis of the numerical results we have obtained, beginning 
with results for the structure of the chains, proceeding to 
the dynamical matrix, and concluding with an analysis 
of the damping of modes in the systems. Section IV gives 
our final conclusions, while certain technical details are 
presented in three appendices. 



II. METHOD 

As will become evident in the forthcoming discussion 
it is advantageous to use two different ways to label the 
atoms forming the junction; these two schemes are illus- 
trated in Fig. [TJ The first scheme (Fig. [T]top panel) is 
based on the cross-sectional area and collects all atoms 
with equilibrium positions on the one-dimensional line 
joining the two surfaces into a "Chain", and calls the 
remaining atoms between the Chain and the substrate 
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FIG. 2: Parameters used for calculating the dynamical ma- 
trix. Only nearest-neighbor coupling is shown. The coupling 
elements labeled 'EM' are found by the empirical model and 
the coupling elements labeled 'DFT' by DFT. On-site ele- 
ments are determined from the coupling elements (see Sec. [A} 
and are not shown in the figure. 

the "Base". The second scheme (Fig. [IJ bottom) dis- 
tinguishes between a "Pyramid" and a "Central Chain"; 
this is chosen because the last atom of the Chain has 
bonds to four or five atoms making this atom very dif- 
ferent from the Central Chain atoms that only have two 
bonds per atom. 

A quantity of central importance to all our analysis 
is the mass-scaled dynamical matrix, K, which we here 
define as including ft, 
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where E is the total energy of the system, Ui is the coor- 
dinate corresponding to the i'th degree of translational 
freedom for the atoms of the system, mi is the mass of 
the atom that the i'th degree freedom belongs to. K gov- 
erns the the evolution of the vibrational system within 
the harmonic approximation. In the Fourier domain the 
Newton equation of motion reads 

Kit a = 4«A i (2) 

where A denotes a mode of oscillation in the system and 
e\ is the corresponding quantization energy. 

The evaluation of K proceeds as follows. Finite differ- 
ence DFT calculations (for details of our implementation, 
see App. [A| were used for the Chain, the Base and the 
coupling between the surface and the Base while for the 
surfaces we used an empirical model due to Treglia and 
Desjonqueres^. Figure [2] illustrates the domains for the 
two different methods. The position of the interface be- 
tween the region treated by DFT and the region treated 
by the empirical model is a parameter that can be varied, 
and the dependence on the final results of the choice of 
this parameter is analyzed in App. IB 21 

The empirical model can be used to describe the on- 
site and coupling elements of atoms in a crystal structure. 
The model uses the bulk modulus of gold to fit the varia- 
tion in the force constant with distance between nearest 
and next-nearest neighbors. Even though the empirical 
model is fitted to the bulk modulus, which is a low fre- 
quency property, it still accurately predicts the cutoff of 
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the bulk band. The positions of the neighboring atoms 
can only have small deviations from perfect crystal po- 
sitions (eg. bulk, surface and ad-atoms). Note that this 
model is general enough to give different coupling ele- 
ments between surface atoms and bulk atoms. The model 
also distinguishes between the coupling between surface 
atoms with or without extra atoms added to the surface. 

The DFT calculations were done with the SIESTA 
code, using the Perdew-Burke-Enzerhof version of 
the GGA exchange-correlation potential with stan- 
dard norm-conserving Troullier-Martins pseuodopoten- 
tials. We used an SZP basis set with a confining energy of 
0.01 Ry. A mesh cutoff of 150 Ry was used. Relaxation 
was done with a force tolerance of 0.002 eV/A. These 
values were found to have converged for the same type of 
system by Frederiksen et ah 31 ^ 32 . The experimental fee 
bulk lattice constant of 4.08 A was used. Only the device 
region was relaxed (defined in Fig. \3§ . 

For the (111) orientations a 4x4 atom surface unit and 
a 2x3 fc-point sampling was used, while for the (100) 
orientations a 3x3 atom surface unit cell and a 3x3 k- 
point sampling was employed. This ensured a similar and 
sufficient fc-point density for both kinds of surfaces (see 

App.EU. 



A. Green's function for a perturbation on the 
surface 

All properties of interest in the present context can be 
derived from the (retarded) Green's function D, defined 

by 

[(e|i)j) 2 I-K]D(e) = lEMD , (3) 

where n = + and we defined the inverse of the Green's 
function by M = D™ 1 . Specifically, we shall need the 
Green's function projected onto the region close to the 
Chain. Our procedure is based on a method due to Mingo 
et alM- which has previously been tested in an investiga- 
tion of finite Si nanowires between Si surfaces. We define 
~X.yz as the the block of the matrix X, where the indices 
run over the degrees of freedom in regions Y, Z, respec- 
tively, where Y, Z = {1, 2, A, D, L, R], as defined either 
in Fig. [l]or Fig.H 

First, let us start with two perfect surfaces. We then 
add the atoms that connect these surfaces (the Base and 
the Chain) . Within a certain range from the added atoms 
the on-site and coupling elements of K will be different 
from the values for the perfect surface. Together, the 
added atoms and the perturbed atoms define the device 
region D (Fig. O bottom). The coupling between the 
device region and the rest of the surface (L, R for the 
left and right leads, respectively) is assumed to be un- 
perturbed. 

In order to compute the Green's function projected on 
the device region, T)dd(^), we first consider this matrix 



Unperturbed 




FIG. 3: Adding atoms to two surfaces. [Top] The forces 
between surface atoms within next-nearest neighbor dis- 
tance^. 08A) of the added atoms are perturbed by the pres- 
ence of the added atoms. [Bottom] The device region is where 
the coupling between the atoms is different from the values 
for the two unperturbed surfaces. The coupling between the 
device region and the leads is considered to be unperturbed. 



representation of Eq. ^ 49 : 

( M DD M Da \ ( T) DD V Da \ = ( I DD Da \ 
\ M aD M aa J \ T> aD T> aa J \ aD l aa J ' 

(4) 

Here the index a = (L, R), i.e., the left and right unper- 
turbed surface, while D = {1,^4,2}. Using straightfor- 
ward matrix manipulations one finds 

Hdd = [M DD -M D a(M aa )- 1 M aD ]~ 1 

= [Mdd-Udd]- 1 , (5) 

which defines the self-energy Hdd = 
MD Q (M aQ )~ 1 M Q £). Since the added atoms do 
not couple to the unperturbed surfaces, and the per- 
turbed region 1 couples only to the right unperturbed 
surface while the perturbed region 2 only couples to the 
left unperturbed surface, the self-energy Hdd has the 
matrix structure 

/M 1L (M LL )- 1 M L1 \ 

U DD = 

\ M 2R {M RR )- 1 M R2 J 

(6) 

This object can be evaluated as follows. First, in the 
limit of large regions 1 and 2, the coupling elements Mli 
and must approach those of the unperturbed sur- 

face, Mjrj and M fl2) respectively. In what follows, we 
shall make the approximation that the regions 1 and 2 
are chosen so, that this condition is satisfied sufficiently 
accurately. Second, we note that the matrix M QQ is in- 
distinguishable from the matrix Mf Q , as long as the in- 
volved atoms are outside the perturbed regions 1 or 2. 
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Therefore, we can write 

M 1L (M LL )- 1 M L1 ~ Mf L (Mf L )- 1 Mf 1 =nf x 
M 2R (M RR )- 1 M R2 ~ Mf fl (M| i? )- 1 M| 2 = n 



22; 

(7) 

where the accuracy increases with increasing size of re- 
gions 1 and 2. On the other hand, using the definition of 
the self-energy, we can write 
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n 



22 



Mf 2 



(D 



22 J 



(8) 



where D 



s 



1,2 is the projection of the unperturbed 



Green's functions onto the atoms in regions 1,2, respec- 
tively. This object is evaluated by exploiting the period- 
icity in the ideal surface plane. The Fourier transform 
of M s in the parallel directions has a tridiagonal block 
structure and we can solve for its inverse very effectively 
using recursive techniques (see e.g. Sancho et alM). Of 
course we still have to evaluate the Fourier transform for 
a large number of Appoints . The density of fc-points as 
well as the size of the infinitesimal i] are convergence pa- 
rameters which determine the accuracy and cost of the 
computation. An analysis of the choice of these parame- 
ters is given in App. IB 21 

To sum up, the calculation is preformed in the follow- 
ing steps: (i) Start with perfect leads and specify the de- 
vice in between them, (ii) The atoms in the leads where 
K is perturbed by the presence of the device are identi- 
fied, (iii) The unperturbed surface Green's function D s 
is found via /c-point sampling and then used to construct 
the self-energy, Eqs.((7H8]). (iv) The perturbed Green's 
function is then found using this self-energy via Eqs.(H- 
i). 



B. Modes and life-times 

For any finite system the eigenvalues e\ and thereby 
also the density of states are found straightforwardly. For 
infinite systems we use that each eigenvector, u\, with 
the corresponding eigenvalue, e\ gives a contribution to 
the imaginary part of the Green's function in the t\ ^> r\ 
limit 
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2e A ( e - £A )2 +7? 2- 
This expression results in the following density of states 



n(e) = lim ImD(e) 

7T jy— >o+ 



(9) 



The broadened vibrational modes of the device region 
can each be associated with a finite life-time. To do this 
we need to have a definition of an approximate vibra- 
tional mode of the central part of the system that evolves 



into an eigenmode of K when the coupling to the leads 
tends to zero. We define 'modes' as the vectors that for 
some energy, e*, correspond to a zero eigenvalue mode of 
ReB DD (e*) (see Sec. O for details). 

We also need to define a few characteristics of a mode. 
The Green's function projected onto a mode can be ap- 
proximated by a broadened free phonon propagator with 
constants e\ and 7a in a neighborhood of the mode peak 
energy, 



u[D DD (e)u x 
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( £ A+7A)+« 2e 7A 



The time-dependent version of the Green's function is an 
exponentially damped sinusoidal oscillation with damp- 
ing rate of mean life-time, T\ = and Q-factor, 
Qx = Comparing the broadened phonon propaga- 

tor to Eq. ([5]) we see that u A ImII(e)?iA = —2ej\, leading 
to 



7A 



4lmII(e*)uA 

2e* 



where e* is the the mode peak energy. 

This calculation of 7a, Qx and t\ only strictly makes 
sense for peaks with a Lorentzian line shape. This re- 
quires that MjImDfl£)(e)iiA is approximately constant 
across the peak which is the case for modes with small 
broadening and large life-time. Nevertheless, we will also 
use these definitions for the delocalized modes since the 
calculated values are still a measure of interaction with 
the leads. 

We also define a measure of spatial localization, s\, 



sa 



E.sclMxl 2 N D -Nc 

T,xeD\c\( U >dx\ 



Nr 



where Nd and Nc are the number of atoms in the device 
and central chain region respectively and D\C means 
Device region except the Central Chain (the perturbed 
reigion). This quantity is useful to pick out modes with 
a large amplitude in the Central Chain region only. We 
have that s\ = 1 signifies equal amplitude in C and con- 
necting atoms, while the limit s\ — ► oo (s\ — > 0) signifies 
a mode which is completely residing inside(outside) the 
Chain. 

It should be stressed that the mode properties calcu- 
lated in this way only refer to the harmonic damping 
by the leads and that other sources of damping are not 
included such as electron-hole pair creation and anhar- 
monicity. The damping due to electron-hole pair cre- 
ation, obtained by an ab-initio calculation on a selection 
of gold chains, is about 50 — 80 /xeV for the vibrational 
mode with the strongest coupling to electron a 30 ' 31 . This 
type of damping is less dependent on strain in gold chains 
due to stable electronic structure as evidenced by the 
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robust electronic conductance of one conductance quan- 
tum. The harmonic damping due to the leads is typically 
higher than this, but as we shall see it can actually drop 
well below this value and thus be less than the electron- 
hole pair damping. 

In case of an applied bias the high frequency modes 
may be excited to a high occupation. The creation of 
vibrational quanta is roughly proportional to eV — huj\, 
while the damping mechanisms are not expected to have 
a strong dependence of the bias. Therefore, as the bias 
is increased beyond the phonon energy threshold, the 
mode occupation will rise and anharmonic interactions 
may become increasingly important even for low tem- 
peratures. Mingo^ has studied anharmonic effects on 
heat-conduction in a model atomic contact, and more 
recently Wang and co-workers et alM- has used ab-initio 
calculations to access the effect of anhamonicity on heat- 
conduction in carbon-based systems. Anharmonic effects 
are, however, outside the scope of the present work. 



III. RESULTS 



A. Geometrical Structure and the Dynamical 
Matrix 



In this subsection we investigate the geometrical struc- 
ture of the chains and the behavior of the dynamical ma- 
trix. For each type of calculation (identified by the num- 
ber of atoms in the Chain, the surface orientation and 
the type of Base) a range of calculations were set up with 
the two surfaces at different separations, Li (i = 1,2...), 
with the separations incremented in equally spaced steps. 
Trial-and-error was used to determine suitable step sizes 
for the different types of calculations. 

To be able to compare chains of different lengths and 
between different surfaces we define the average bond 
length, B = (bj), as the average length between neigh- 
boring atoms within the Chain, where j runs over the 
number of bonds in the Chain (see Fig. [4]). B is useful 
because it is closely related to the experimentally mea- 
surable forced on the Chain and can be found without 
interpolation. The close relationship between B and the 
force is demonstrated in Fig. [5] where the force is calcu- 
lated as the slope of a least-squares fit of 

[{Ei-l,Li-i), (Ei, Li), (JEj+i,£j+i)] 

where E is the total energy. We note that the force vs. 
B curves to a good approximation follows a straight line 
with a slope of k = 2.5 eV/A , which can be interpreted 
as the spring constant of the bonds in the Chain. In 
addition to B we also define the average bond angle T = 
(0 S ). 

The behavior of the systems with respect to B is rel- 
atively simple. As the systems are strained it is mostly 
the bonds in the chain that are enlongated. Finally the 
central bond(s) become so weak that they break. At 




FIG. 4: Distances and angles used to define the average bond 
length, B = {bj}, and the average bond angle, T = (9j), 
respectively 
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FIG. 5: (color online) Force as function of average bond 
length, T = (0j). 



low B we see from Fig. [6] that the longer chains adopt 
a zig-zag confirmation at low average bond length. The 
3- and 4-atom chains, however, remain linear within the 
investigated range. Furthermore, the longer chains have 
a similar variation in the average bond angle. 

These preliminary observations are in agreement with 
previous theoretical studies by Frederiksen et alM- and 
Sanchez-Portal et aZ.— . We recount these observations 
because we find that using B as a parameter provides a 
helpful way to compare chain of different lengthts and be- 
cause the calculations in this paper are the most accurate 
to date££. 

To shed light on the effect of straining the chains, we 
next investigate the energies that are related to differ- 
ent types of movement by analyzing the eigenmodes and 
eigenvalues of selected blocks of K. Especially, we can 
consider the local motion of individual atoms or groups 
of atoms, freezing all other degrees of freedom, by pick- 
ing the corresponding parts of K. For a single atom this 
amounts to the on-site 3x3 blocks. The square root 
of the positive eigenvalues of the reduced matrix, which 
we call local energies, give the approximate energy of a 
solution to the full K that has a large overlap with the 
corresponding eigenmode, if the coupling to the rest of 
the dynamical matrix is low. The negative eigenvalues 
of a block are ignored since they correspond to motion 
that is only stabilized by degrees of freedom outside the 
block. 

The behavior of the dynamical matrix in terms of lo- 
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FIG. 6: (color online) Average bond angle as a function of 
the average bond length. The long chains adopt a zig-zag 
structure at low B while the short chains remain linear. 



20 



> 

gli 



- 

"Bo 

I 1 

c 

w 

CJ 

o 
- 



B=2.65 A 




B=2.88 A 




4 Atoms 
(100) 



I 11 



surf, base 



1 



2 3 4 base surf. 



FIG. 7: (color online) Local energies of a 4 atom chain be- 
tween two (lOO)-surfaces at different strains. The largest 
eigenvalues are connected by a line to guide the eye. 



FIG. 8: Local energies for selected blocks of K (Pyra- 
mid/Central Chain) plotted vs. the average bond length in 
the 4 atom chain. Since the central chain in this case consists 
of 2 atoms we can classify the eigenvectors as LO: longitudi- 
nal optical, LA: longitudinal acoustic, TO: transverse optical 
(doubly degenerate) or TA: transverse acoustic (doubly de- 
generate). 



B. Mode life-times and Q-factors 

We next investigate the modes of the finite chain sys- 
tems. An example is given in Fig. [9] which depicts the 
projected DOS for a chain with 4 atoms at an intermedi- 
ate strain. Notice the large variation in the width of the 
peaks. The peaks with a low width correspond to modes 
that have the largest amplitude in the Chain, while the 
peaks with a large width correspond to modes with large 
amplitude on the Base and Surface. Since this type of 
system has no natural boundary between 'device' and 
'leads' we will have large variation in the harmonic damp- 
ing no matter where we define such a boundary. 



cal energies, is relatively straightforward, as illustrated in 
Fig. \7\ When the bonds are strained they are also weak- 
ened. In the Central Chain the local energies are quickly 
reduced with increased strain (w 65% decrease) while 
the dynamical matrix of the surfaces is hardly affected. 
The Base and the first atom of the Chain fall in between 
these two extremes with a 20% and 40% decrease, respec- 
tively. The middle bonds in the Central Chain are the 
ones that are strained and weakened the most when the 
surfaces are moved apart. It is also where the chain is 
expected to break^. Most interestingly, we note that at 
least one jump in the on-site local energies occur when 
moving from the Surface to the Central Chain. 

In Fig. [8] we see how motion parallel to the chain is at 
higher energies than perpendicular motion, and that the 
LO type motion of the ABL modes has the highest en- 
ergy. We also see that the local energies of the ABL/LO 
type motion moves past the local energies of the Pyra- 
mid as the strain is increased. In this way the ABL/LO 
modes can in some sense act as a probe of the contacts. 
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FIG. 9: Projected DOS onto a representative selection of the 
vibrational modes of the device region that has a large overlap 
with the added (Chain+Base) region. 

In Fig. [TO] we present the Q-factor, spatial localization 
and peak energy of all modes for chains with 3-7 atoms 
between (100) surfaces. These are the main result of 
this article. Table [TJ shows the same information in an 
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alternative form. We now proceed to an analysis of these 
results. 

The ABL modes are of special interest. These 
modes have been identified by previous theoretical 
and experimental studies as the primary scatterers of 
electronfl!2i2iL2!L2JL22isyyi. The ABL modes are easily 
identified in Fig. [10] since they have the highest energy 
of the modes that are spatially localized to the central 
chain (black or dark gray on the figure). Modes corre- 
sponding to transverse motion of the central chain are 
also clearly visible. These modes are energetically and 
spatially localized, but are of limited interest because of 
a low electron-phonon coupling. 

Certain ABL modes are very long-lived. At low strains, 
ABL modes lie outside the bulk (and surface) band and 
have, in our harmonic approximation, an infinite Q- 
factor. In reality the Q-factor will be limited by electron- 
phonon and anharmonic interactions. At higher strain 
the ABL modes move inside the bulk band and one ob- 
serves a great variation in the corresponding Q-factors. 
When the peak energy lies inside the bulk band there ex- 
ists modes in the bulk with the same energy and it will 
mostly be the structure of the connection between the 
bulk crystal and the chain that determines the width of 
the peak. 

The long chains tend to have longer lived ABL/LO 
modes due to the larger ratio between the size of the 
Central Chain and the size of its boundary. The 7-atom 
chain is especially interesting since it has an ABL/LO 
type mode with a damping of 5 meV at one strain, while 
at another strain the ABL/LO mode has a damping of 
300 meV, i.e., more than one order-of-magnitude varia- 
tion in the harmonic damping of the primary scatterer 
of electrons due to only a 0.03 A change in the average 
bond length! 

The largest damping of an ABL-mode for these sys- 
tems is 7a ~ 1 meV, which is still significantly lower 
than the w 20 meV band width. This can be attributed 
to fact, noted above, that there always exists a large mis- 
match in local energies moving from the central part of 
the chain to the rest of the system (see Fig. [7j) . 

Previous studies by Frederiksen et al£^ obtained 
a rough estimate for the variation of the non- 
electronic (harmonic and anharmonic) damping of 5- 
50 /xeV for the longer chains by fitting the experimen- 
tal IETS signals of Agrait et alr^ to a model calcula- 
tion. The estimated peak energies lie well within the 
bulk band for all the recorded signals. The reason the 
non-electronic damping rate can be extracted is because 
the excitation of vibrations and damping of vibrations 
through electron-hole creation are both proportional to 
the strength of the electron-phonon coupling. This means 
that the step in the experimental conductance, when the 
bias reaches the phonon energy, can be used to estimate 
strength of the electron-phonon interaction and thereby 
the electron-hole pair damping. The slope in the con- 
ductance beyond this step can then be used to extract 
the total damping. By subtracting the electron-hole pair 



damping from the total damping we get an estimate of 
the sum of the sum of harmonic and anharmonic contri- 
butions to the damping. 

The estimate in Ref. [3l| agrees well with our lowest 
damping of 5 /J,eV. The highest damping we have found 
was w 400 /zeV found for the 6 atom chain which is an 
order of magnitude larger than the upper limit of Ref. 
[3lj . We believe that this discrepancy can be largely at- 
tributed to the difficulty in extracting the necessary pa- 
rameters from experiments when the harmonic damping 
is large. Furthermore, for the 6-7 atom chains we ob- 
serve that the high damping occurs at low strain, where 
the electron-phonon coupling is weakii. 



Chain 


Qx 


7A(^eV) 


ta(ps) 


3(100) 


3-7 


800-1200 


0.5-0.8 


4(100) 


5-30 


100-900 


0.7-7 


5(100) 


10-40 


200-500 


1.3-3 


6(100) 


15-80 


90-400 


1.6-7 


7(100) 


40-1500 


5-300 


2-130 


5(111) (symmetric) 


15-80 


40-400 


1.6-16 


5(111) (asymmetric) 


10-100 


40-800 


0.8-16 



TABLE I: The variation of the Qx, 7a and r A of the ABL/LO- 
modes. For Chains with 3-7 between (100) surfaces and for 
Chains with 5 atoms between (111) surfaces but with slightly 
different Bases. The strains where the peak energy of the 
ABL/LO-mode falls close to or outside the bulk band edge 
have been disregarded. 

There are two main differences between the (100) and 
the (111) systems. The first difference is that the (111)- 
systems have ABL/LO-modes that are more long-lived 
compared to the (lOO)-systems (see TableUand Fig. fTT]) . 
The second difference is the behavior of the localized 
modes close to the band edge (See Fig. [TTJ. The modes 
with energies outside the bulk band in the (111) systems 
are less spatially localized compared to the (100) case. 
At low strain, the (lll)-chain have ABL/LO-modes ex- 
tending further into the Base and Surface than the (100)- 
chain. 

There are certain general features of how the damping 
evolves with strain that are easily understood. Modes 
with peak energies in the range 16—19 meV in generel 
have a very high damping while those in the range 
14—16 meV have very low damping. This correlates 
well with the bulk DOS for gold (see e.g. Ref. [1§|). 
The optical peak in the bulk DOS corresponds to strong 
damping while the gap between optical and acoustical 
modes correspond the the range of low damping. 

To sum up, localized modes occur at low strain where 
the bonds in the chain are very strong, and give rise to 
frequencies close to or outside the bulk band edge. Inside 
the bulk band strong localization is still possible for the 
long chains, especially the 7-atom chain. This requires, 
however, that the coupling between the Central Chain 
and the surface is weak at the typical frequency of the 
ABL/LO mode due to the structure of the connection. 



8 



20 - 



3 Atoms 
Bulk Rand Edge ( 10 °) 



Q=7 




i 



Q=3 



4 Atoms 
(100) 



■ • 

t ( 



T f i i i 



i i t i 
itii 
> i i i 




Q=30 



1 1 1 1 1 1 
.., 5 Atoms 

ft, \ (100) 


i 


r T! I ! 

* t • 

; t 

i i It 

i III - 

• 

i i.i i 


- 
- 

Q=10 

■ 

» 

* 



i i i i i 1 

6 Atoms 

a Q= 80 (ioo) 


* \ 

? 

■ 

- 'til 
111! 

' " » r 

: : : i 

: ; j | ; 

• • i • 

• * • • 


I i 

♦ i 

' i ! 

k * i 

• • 

• 


- 

Q=15 



1 








i 

7 Atoms 


F • 


Q= 


=40 




(100) 


: 

- 


t 


T 


i i 

i I 




- 

: 
: 


• 

: 

; 

: 


• 

♦ 

! 


• 

i * 


Q=1500 

• 

t 


t 
i 


t 

i 
i 


t 

? 


t 

t ' 


t 



2.65 2.70 2.75 2.80 2.85 2.90 

B(A) 

FIG. 10: The vibrational modes for chains with 3-7 atoms be- 
tween two 100-surfaces. The center of the disks are positioned 
at the peak of the projection of vibrational DOS on the mode 
in question. The area of a disk is proportional to the Q\, but 
is limited to what corresponds to a Q-factor of 250. he gray 
level, that ranges from light gray to black in 4 steps signifies 
that s A G [0, 2 [(light gray), [2, 4[, [4, 6[, [6, 8[ or [8, oo [(black). 
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FIG. 11: The vibrational modes for 5 atom chains between 
two (111) surfaces. [Top] Symmetric pyramids. [Bottom] 
Asymmetric pyramids (one atom added to one of the pyra- 
mids). The area of a disk is proportional to the Q\, but is 
limited to what corresponds to a Q-factor of 250. The gray 
level, that ranges from light gray to black in 4 steps signifies 
that s A € [0, 2[(light gray), [2, 4[, [4, 6[, [6, 8[ or [8, oo [(black). 



The behavior depends strongly on the detailed structure 
of the base and the state of strain, but some general 
features can be related to the the bulk DOS. 



IV. CONCLUSION AND DISCUSSION 

We have presented a study of the harmonic damp- 
ing of vibrational modes in gold chains using a method 
that uses ab-initio parameters for the chains and em- 
pirical parameters for the leads. We have focused on the 
ABL/LO modes that interact strongly with electrons and 
are thereby experimentally accessible through IV spec- 
troscopy. We provide an estimates for the damping of 
ABL/LO-modes from ab-initio calculations as a function 
of strain for a wide range of gold chain systems. The cal- 
culations of the ABL-phonon damping rates agree well 
with earlier estimates, found by fitting a model to exper- 
imental inelastic signal a n i 32 . 

We have found the that the values of the harmonic 
damping for the ABL modes can vary by over an order 
of magnitude with strain. Even with small variations in 
the strain, the harmonic damping can exhibit this strong 
variation. This extreme sensitivity may explain the large 
variations seen experimentally in different chains. 

The range of the harmonic damping also depends 
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strongly on the number of atoms in the chain since we 
see a clear increase in localization going from a 6- to a 
7-atom chain. The chain with 7 atoms really stands out, 
since it, in addition to having very localized modes in 
generel, it also has the greatest variation in harmonic 
damping. This strong variation in the harmonic damp- 
ing of the ABL/LO-modes, that depends on the details 
of the structure, suggest that accurate atomistic calcula- 
tions of the vibrational structure is necessary to predict 
the inelastic signal. 

All types of chains were found to have ABL/LO-modes 
tha lie outside the bulk phonon band at low strain. These 
modes are expected to have very long life-times since 
the harmonic damping is zero. Signatures of the rather 
abrupt change in the damping of the ABL/LO-modes 
when strained have not been discussed in experimental 
literature so far. We believe this is due to the com- 
mon experimental techniques for producing these chains 
heavily favor strained chains. The ABL/LO-mode life- 
time may be set by the coupling to the electronic system 
(electron-hole pair damping). Indeed, even inside the 
bulk band the electron-hole pair damping can be of the 
same order as the harmonic damping. For example, a 
7a eh ~ 50 — 80 fieV was found for a 4-atom2& and a 7- 
atom^I chain, which we can compare with 100 — 900 /ieV 
and 5 — 300/ieV found above for the harmonic vibrational 
damping. Thus the damping can in certain cases be dom- 
inated by the electron-hole pair damping for frequencies 
even inside the bulk band. 



Finally we find a difference in the the damping of 
ABL/LO modes in chains between (100)- and (111)- 
surfaces. For the investigated 5-atom chains there is both 
a marked difference in the strength of damping and in the 
variation of the damping with strain. It might be possible 
to distinguish between (100) and (111) pyramids exper- 
imentally due to this difference. The ABL/LO modes 
will have strong coupling to the bulk at certain ener- 
gies, characteristic of the pyramid type. This in turn, 
results in broadening/splitting of the modes depending 
on whether the characteristic energies are inside or out- 
side the bulk band. This broadening/splitting would be 
detectable in the TV-curve since it is related to the char- 
acteristics of the conductance step at the peak energy 
of the vibrational mode. Finding IV-curves at different 
strains could thereby serve as a fingerprint of the specific 
way the chain is connected to the surroundings. Hihath 
et al. have demonstrated that such measurements are in- 
deed possible on a single-molecule contact. 

The techniques used in this paper can be combined 
with electronic transport calculations to predict the in- 
elastic signal in the IV characteristic of a system. This 
will be done in future work, where we will also elimi- 
nate the use of the empirical model for the leads and use 
ab-initio parameters for the entire system. 
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APPENDIX A: CONSTRUCTING THE 
DYNAMICAL MATRIX 

In this subsection the details of how we constructed 
the dynamical matrix are presented. The dynamical ma- 
trix must be symmetric and obey momentum conserva- 
tion. Momentum conservation, in this context, means 
that when an atom is displaced the force on the dis- 
placed atom equals minus the total force on all other 
atoms. We ensure momentum conservation by setting 
the on-site 3x3 matrix to minus the sum of the force 
constant coupling matrices to all the other atoms. This 
method for regularizing the dynamical matrix was previ- 
ously used by Frederiksen et aL—, and generally improves 
on the errors introduced in the total energy when displac- 
ing atoms relative to the underlying computational grid 
(the DFT egg-box effect). We calculate off-diagonal cou- 
pling part of the force constant matrix was calculated 
with a finite difference scheme using a displacement, Z, 
of 0.02 A in the x,y and z directions for all atoms in the 
Chain and Base. 

To improve the accuracy further, the forces were cal- 
culated for both positive and negative displacement. If 
i and j are degrees of freedom situated inside the DFT 
region we therefore perform 4 independent calculations 
of Ky = Kjj, since K is a symmetric matrix. In the 
end we use the average of the force constant from these 
4 calculations 



F, 



Z 



Z 



%)/4, 



where e.g. Fy + denotes the force on i due to a positive 
displacement of j. If i is inside the DFT region and j is 
not, the coupling is calculated as an average of 2 force 
constants 



K, 



,F 



^/miirij 



F 



ji- 



0/2. 



If an atom was close to a periodic image of another 
atom (less than half the unit cell length in any direction) 
the force between these atoms was set to zero to avoid 
artifacts of the periodic calculational setup. The empiri- 
cal model was used to calculate the coupling between the 
surface atoms. After all coupling elements were found 
the on-site elements were calculated for the system as a 
whole. 
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APPENDIX B: CONVERGENCE 
1. Convergence Parameters 

In the calculations there are several convergence pa- 
rameters, and here we provide an overview. 

There are three important length-scales in the calcu- 
lations: L\,L<x and L3. We assume that when two atoms 
are further apart than L±, the coupling elements between 
them vanishes. L 2 is the correlation length for proper- 
ties that do not have an energy-dependence, like forces, 
equilibrium positions and total energies, while L3 is the 
assumed correlation length for properties that do have 
an energy-dependence, like the surface Green's function, 
vibrational DOS etc. L3 always needs to be larger than 
L2, L3 > L/2, but the specific size needed depends on 
the required energy resolution. L 2 determines the fc- 
point sampling used in the DFT-calculations and L3 the 
fc-point sampling used in the calculation of the surface 
Green's function. In each case the number of fc-points 
used one direction is chosen to be the smallest integer, 
i, such that i > ^) where a is the size of the calcula- 
tional cell in that direction. The DFT fc-point sampling 
used is dense enough to ensure that L 2 > 23 A for all 
calculations. 

In the calculation of the Green's functions we intro- 
duced a finite artificial broadening. This broadening, 77, 
was divided into a small broadening of the device region, 
77c, and a large broadening for the leads, t]l- The rea- 
soning behind this is that the density of states is much 
more smooth in the bulk-like regions far away from the 
chain. A large tjl has the advantage that it reduces the 
need for fc-point sampling drastically. Without a small 
r\c we would not be able to discover very sharp peaks 
in the DOS. To reliably find the modes of the system it 
is also important that the energy spacing, AE is on the 
same level or smaller than r\c- 

The artificial broadening limits how large life-times we 
can resolve. This is why we in the following write the 
upper limit to the life-time introduced by the artificial 
broadening. 

A final convergence parameter is the position of the in- 
terface between DFT and empirical model parameters for 
the dynamical matrix. This is a very important param- 
eter since the error introduced by having this interface 
relatively close to the chain is what limits the precision 
of the calculations. 



2. Test of Convergence 

Next we present the tests that have been carried out to 
ensure that the calculations in this article are sufficiently 
converged. The convergence for the SIESTA basis set 
and the size of the finite displacement used in the finite 
difference calculations was already tested for the same 
type of systems by Frederiksen et aZ*2I. 
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FIG. 12: Modes of the device region with a different DFT/EM 
interface. The label designate the region treated with DFT, 
where the 'Added' region is the one used in the main part 
of the calculations and 'AH' is fully ab-initio. See Fig. 1111 to 
see what the color and size signify. On this plot modes with 
S\ € [0, 2[ are suppressed. 



So here we first examine the convergence of the DFT 
calculations of the dynamical matrix. A calculation for 
a 4 atom chain between (lOO)-surface was done with im- 
proved values for the important DFT convergence param- 
eters. The mesh cutoff was increased from 150 to 200 Ry, 
and the fc-point sampling was increased from 2 x 3 to 
3x4. For this change in parameter we obtained a max- 
imal difference of 0.2 meV, when comparing the square 
root of the sorted array of eigenvalues of the dynamical 
matrix. This is a negligible size since the average value 
of the eigenvalues is about 10 meV. The fc-point sampling 
in the DFT calculations proved crucial for the structure 
of the strained systems, since gamma-point calculations 
resulted in different structures (different bonds weakened 
at high strain) with very large life-times. 

The perturbation length used in our calculations was 
L\ = 4.08A, which is the same as next-nearest neighbor 
distance. The magnitude of any next-nearest-neighbor 

coupling matrix, defined as |X| = ^Jj2ij ~%fj was never 

larger than 15% compared to the magnitude of any 



nearest-neighbors coupling matrix. The error introduced 
by this truncation is smaller than the one introduced by 
using the empirical model for the dynamical matrix. 

For the calculation of the DOS we gradually im- 
proved L3, rji and rjc and found that the DOS was 
converged using AE = 10 /LteV, r\L = 100 /jeV(7 ps), 
?yc = 10 /jeV(70 ps) and L 2 = 200 A(68x68 fc-points) 
except in one calculation for the 7 atom chain we needed 
the life-time of one very sharp peak. This required a bet- 
ter resolution using AE = 1 /LteV, t\l = 10 /_teV(70 ps), 
ric = 1 MeV(700 ps) and L 2 = 400 A (136x136 fc-points). 

Finally, we have considered how much the interface be- 
tween the DFT parameters and the empirical parameters 
affect our results. In Fig. [12] we show a study where we 
vary the position of this interface. We find that our cal- 
culation of the Q-factor and the spacial localization is 
converged to about the first significant digit for modes 
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FIG. 13: Example ReD and n(e), dashed and solid line re- 
spectively, for a Green's function with two poles at 1 and 2 
with a 0.1 broadening. We see that the values where the real 
part is zero only correspond to peaks in the density if the 
slope positive. 



quirement, for the defition of modes in the case of the 
open system, is that these modes become the modes of 
the isolated system in the limit of zero coupling between 
the device region and the leads. 

The following definition fulfills this condition. A 
'mode' is defined as a (complex) eigenvector u\ of 
M Dfl (f*) (and D DD {e*)) that fulfills, 



Re{ulD DD (e*)u x } = , (CI) 



and 



—Rc{uln DD {e)u x }\ e=e ,>0 (C2) 



that are spatially localized to the Central Chain. We 
judge that this is what mainly sets the limit of accuracy 
of in our calculations. 



APPENDIX C: DEFINITION OF THE MODES 
FOR AN OPEN SYSTEM 

The starting point is the modes of a closed system, 
namely, the eigenmodes of K. The most important re- 



fer some energy, e*, corresponding to a peak in DOS. 

An illustration of these two conditions is given in 
Fig. [13l In practice, the modes are found from the num- 
ber of positive eigenvalues of Ddd evaluated at each 
point of our energy-grid. If this number increases be- 
tween two successive energies, e and e + Ae, the eigen- 
modes at these two energies are matched up. The eigen- 
mode corresponding to the eigenvalue that changes sign 
is then identified as a mode of the open system. 
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